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The presence of spiral structure in isolated galaxies is a problem that has only been partially 
explained by theoretical models. Because the rate and pattern of star formation in the disk must 
depend only on mechanisms internal to the disk, we may think of the spiral galaxy as a self-regulated 
system far from equilibrium. This paper uses this idea to look at a reaction-diffusion model for the 
formation of spiral structures in certain types of galaxies. In numerical runs of the model, spiral 
structure forms and persists over several revolutions of the disk, but eventually dies out. 
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I. INTRODUCTION 



The problem of how spiral structures form in galaxies is one that has often been studied (for a general overview of 



the subject, see [0), and can be divided into two aspects. The first is that of temporary structures most likely caused 
by gravitational perturbations from passing galaxies or an asymmetric halo, or disk material having an initial velocity 
relative to the local standard of rest from formation processes. Here the phenomenon seems to be related to density 

■ waves H , quasi-stable modes in the gravitational potential of the disk which offer a good description of grand design 
spirals, where the arms are well-defined and has a high degree of symmetry. However, this paper is devoted to the 

• second part of the problem, namely, that of spiral structures in "isolated" galaxies, i.e. where we ignore influences 
from outside the disk and consider only internal (and recurring) mechanisms. These galaxies can be represented 
by flocculentj^] spirals, so called because of the fleecy appearance of their many short and asymmetric spiral arms. 
Q^' Because these spirals are seen in blue light, but not red || - suggesting the arms are composed of younger, bluer stars, 
while the older, redder stars are more evenly distributed across the disk - they must be an artifact of star formation 

■ itself, and are not primarily density waves. Our attitude here is to consider structure formed by the processes of star 
formation and neglect gravitational influences. This should apply equally to all spirals, although for grand design 

C$ • spirals, we would expect it to act with density waves arising from other sources. 

Most of the theoretical work done in star formation processes deals with the solar neighborhood (those stars in the 
vicinity of our Sun) as opposed to structures on a galactic scale [BJ . Other models used to study these aspects of star 
' formation and spiral structure, such as those built on propagating star formation 0], either greatly simplify the physics 
involved or need finely tuned parameters to match observations. In this paper, we give a model of spiral structure 
in isolated galaxies based on the idea of star formation as part of a network of self-regulated and autocatalyzed 
reactions f?],||]0. Recently, there has been a great deal of both theoretical and experimental work studying how non- 
equilibrium systems in chemistry and biology produce patterns in both time and space. This has included looking at 
both the organic - such as bacterial colonies , the differentiation of cell types [Efl , and the formation of embryonic 



structure in multicellular organisms |10| - as well as the inorganic - the BZ chemical reaction diffusion limited 
aggregation jl^ ], and self-organized critical systems p3[ . There have been many successes in reproducing patterns 
in the laboratory, and these models typically share both partial differential equations and discrete elements such 
as cellular automata. However, most of this work in non-equilibrium systems has been after the main work on the 
formation of structures in astronomy. 



1 The division of spiral galaxies into grand design and flocculent is based on the arm classification scheme of Elmegreen and 
Elmegreen ||. 

2 This work, is different from the work || done in the past, because it uses a modified model. In particular, the models in [^] 
had a subtle flaw in normalizations and were also missing some physically important reaction terms. See for a full discussion 
of this. 
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Thus one might seek to apply this work to patterns seen in spiral structures. Because we are studying the process in 
isolated galaxies, we know that the formation must be caused by events within the disk, rather than by the actions of 
outside players. The isolated galactic spiral is far from equilibrium - there is differentiation of material into stars and 
clouds of gases whose distribution varies over space and time. In addition, star formation happens at a constant rate[], 
as averaged across the disk. This is a clue that the process is regulated by a feedback loop to maintain this constancy 
(for evidence of this mechanism, see fig] ). These characteristics are shared by other types of non-equilibrium systems. 
Below, wc list the predominant features that these kinds of networks of reactions have in common, along with examples 
of the same behavior in galactic disks. 

• Steady state system There is a slow (relative to the dynamical time scales) and steady flow of energy, 
and perhaps matter running through the disk. In spiral galaxies, star formation proceeds at a constant rate, 
averaged over the disk, for time scales on the order of 10 10 years. The fact that this is greater than the time 
scales of the actual star formation process (10 7 years) implies that the slow and steady rate is regulated by 
feedback mechanisms. 

• Non-equilibrium system The steady state is far from thermodynamic equilibrium, and there is a coexistence 
of several species or phases of matter, which exchange matter and energy among themselves through closed cycles. 
The galactic disk is not a uniformly dense clump of material at thermal equilibrium, but instead is divided into 
gases at different temperatures and stars of various mass. These species exchange matter: for example, massive 
stars supernova to form warm gas, which can cool and then condense into new stars. 

• Feedback mechanisms The rates at which material flows around these cycles are governed by feedback loops 
that have arisen during the organization of the system into the steady state. An example of this is suggested 
by Parravano and collaborators fil| which explains how the average pressure in the interstellar medium (ISM) 
is maintained. They argue that there are two phases, the warm gas of the ambient phase, and the cooler gas of 
the condensed phase, with a phase boundary in the pressure-temperature plane. Ultraviolet radiation from the 
supernovae of massive stars heats the gas, which prevents the condensation of newer stars, so the supernova rate 
goes down, allowing new massive stars to form (and so the supernova rate will increase again). This feedback 
mechanism keeps the gas on the phase boundary. 

• Autocatalytic reaction networks Any substances that serve as catalysts or repressors of reactions in the 
network are themselves produced by reactions inside the network. Suppose we look at the condensation of giant 
molecular clouds (GMCs). This is catalyzed by dust grains produced by cool giant stars, shielding the clouds 
and providing sites for molecular binding, and carbon and oxygen, which may cool the clouds by radiation from 
the rotational modes of CO molecules. The condensation is inhibited by ultraviolet radiation from massive stars, 
as described by the Parravano process mentioned above. 

• Separation in space There may be spatial segregation of the different phases or materials in the cycles. 
This occurs when the inhibitory and catalytic influences propagate over different distance scales. At the smallest 
scale, this means the production of certain substances may be subject to refractory periods - once production 
has occurred in a local region, it will not be repeated there for some period of time. For the influences in 
the process of GMC condensation, dust grains, carbon and oxygen propagate only over distances of about 100 
parsecs (how far supernovae and massive stars can spread their products) while UV radiation can travel over 
much of the galactic disk. 

As can be seen, there is evidence that we can think of spiral structure in isolated galaxies to be a product of a self- 
organized, autocatalyzed network of reactions in the star formation process. Given a system with the characteristics 
listed above, there are models which can describe the spatial structure, the most typical of which is the reaction- 



diffusion model p0 17 1, and we develop a model along these lines for spiral structure in isolated galaxies. 



II. THEORETICAL ASPECTS 



3 This is true in observed galaxies up to a factor of two; see Sandage for how the star formation rate in different types of 
spirals changes with time. 



2 



A. One-zone model 



We briefly outline the processes necessary to start star formation. The cold clouds of the GMCs condense out of the 
interstellar medium (ISM), forming distributions of gas and dust that are apparently scale-invariant. As mentioned 
above, the condensation is helped along by the actions of dust, carbon and oxygen, while it can be impeded by UV 
radiation from massive stars. The Parravano process may place some limitations on the amount of condensation. 
Typical time scales for the inhibition are about 10 7 years, the average life span of massive stars, after which the 
supernova (SN) rate and the UV radiation flux will die off. Once the GMCs start to condense, then their cores may 
collapse to form stars. This collapsing is brought on by Shockwaves from supernovae or HII regions (we neglect any 
impact that density waves may have in these collapses), and so have the same length scale as the propagation of dust 
by supernovae, or about 100 parsecs. Once the stars are formed, they can inhibit the infall of gas by the stellar wind 
or UV radiation produced by the star Jig] . These effects occur on short length scales, about the size of one cloud 
complex, and reduce the star forming efficiency of the clouds down to a few percent. 

We must take the processes occurring in the galactic disk and abstract them to produce a viable mathematical 
model. One way we can simplify the system is to take the continuous spectrum of star types and break it up to those 
massive stars that can supernova (and thus provide matter and radiation back into the system) and those lighter stars 
that cannot. We will be neglecting the fact that these lighter stars can return matter to the ISM, so in our model, 
they will simply act as a matter sink. There is also matter exchange between the cold gas of the GMCs and the 
warmer ambient gas, due to heating and cooling. We can summarize the flows of material in Fig. 1, and the various 
material and energy components are given in Table ||. Now we give a brief description of the physical processes we 
include in the model. 

• Cloud destruction cs — » g,cs — > s: Because of the presence of massive stars and their stellar winds and SNe 
produced shocks, there will be a mechanism of cloud distruption. Some of this pressure will trigger star formation 
(cs — > s), but the eventual result will be the disruption of the cloud into warm gas (cs — > g). Typically, the 
efficiency of star formation is around a few percent. 

• Cloud-cloud collisions c 2 — » g, c 2 — > s: Another source of pressure is the collision of clouds, which will have the 
same types of effects as the star-induced cloud destruction mentioned above. Again, the rate of star formation 
is a few percent. 

• Mass infall to the disk 5: Galaxies are believed to be formed from the condensation of matter from a spherical 
halo into a disk, and so there is certainly the possibility that there is a continuing flow of matter. It is believed 
that the rate is enough to replenish the material in the disk in a time span of billions of years (see Section 4.3 

of Larson H). 

• Direct cloud destruction by stars cs — * g: This represents the effects of stars, such as stellar winds, which come 
directly from the massive stars, as opposed to radiation and Shockwaves, which might travel some distance. 
The main physical action behind this term is the ionization and champagne flows of stars formed inside the 
cloud @. 

• UV radiation, Shockwave production s — > r, s — > h: The sources of UV radiation and Shockwaves are from SN 
events from massive stars. These effects are more long range, although shocks will travel only about 100 pc, 
while radiation can traverse the entire galaxy. 

• Damping terms (c + g)r — » r, (c + g)h — * h: Because the energy carried by UV radiation and Shockwaves will be 
dissipated by the interaction with matter - both warm and cold - we include a damping term. 

• Cooling term g 2 /r — > c: UV radiation will act as a thermostat, since warm gas is less likely to cool in an 
environment with a high radiation density. In the results presented here, we consider cooling which is inversely 
proportional to the radiation density. 

• Cloud destruction by UV radiation and Shockwaves cr — ► g,ch — + s: These effects are more long-range than 
the direct cloud destruction used previously. Note that ultraviolet radiation will ionize the clouds into warm 
gas, as does the destruction induced by massive stars, but the pressure due to Shockwaves will initiate new star 
formation. 

Now, we take this model of the material flows in the galactic disk, and write them in a system of equations using 
the law of mass action - that is, the change in time of the output quantity is given by the product of the densities of 
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the inputs, with a constant parameter giving the rate of the reaction. Thus, given the processes outlined above, we 
have 



— = a 2 cr — \5\ch — 7CS — e\c (la) 

at r 

dg otig 2 , 

-7- = ■ h a 2 cr + s + jcs + e 2 c + (lb) 

r 

^ = fcch + e 3 c 2 - s (lc) 

— = 771S - 0i(c + g)r (Id) 
dh , . 

— = 7] 2 s - <p-2[c + g)h (le) 
^=/3 3 ch + e iC ? (If) 

where the parameters of the equations cti, a%, ■ ■ ■ and their ranges are given in Table O. To choose the parameters, 



we use the lifetime of a typical massive star, 10 7 years, as the dimension of time, and the units of mass and energy to 
be those appropriate for each components - for example, we choose the mass unit for warm gas g to be one hydrogen 
atom per cubic centimeter. Then each parameter is the rate of the reaction at the mass and energy densities we 
select, e.g. since the mass flow into the galaxy is estimated to completely replace the current material in 10 10 years, 
we choose the mass inflow constant 5 ~ 0.001. For a fuller discussion of the choices made, see 0. Finally, note that 
conversation of matter implies f3\ — (3 2 + fa and e% = e 2 + e 3 + e 4 . 

As mentioned above, only massive stars are adding material to the ISM, through the mechanism of supernovae - 
we are ignoring the fact that light stars add material via evaporation. There are a few other simplifications which 
have been used to arrive at these equations, such as neglecting the impact of catalyzers to GMC condensation such 
as dust and carbon to avoid parameters that depend on metal concentrations. Also, we noted previously that our 
cooling term is inversely proportional to the radiation density. This gives the simple result that cooling is faster when 
there is less radiation, but this is certainly not the only possibility (although it is the one we have examined more 
closely). Another choice would be to include something like a step function - once the UV radiation is lower than 
a certain amount, the rate of cooling increases greatly, but above the cutoff, it is negligible. This idea of a critical 
density is essentially the process advocated by Parravano and collaborators. However, it can be shown that there 
is little variation in the average values of the components as the functional dependence of the cooling is altered. In 
addition, the use of a step function in the cooling function can lead to oscillatory behavior in the components that is 
almost discontinuous and is therefore undesirable. 



B. Reaction-diffusion model 

The system of equations (|l]) can be useful in understanding such things as the chemical evolution in a galaxy, but 
there already exists a substantial literature on one-zone models in galactic evolution. Our goal here is to work with 
a model, similar to those in chemical and biological non-equilibrium systems, where spatial and temporal patterns 
are generated. In these systems, there are different length scales: L; ong is the scale of the whole disk, the distance 
that UV radiation can travel; Lint is the scale of distances between cloud complexes, and the distance core collapse 
is induced by supernovae; and L s h or t is the scale of a single cloud, the distance new stars evaporate the GMC from 
which they condensed. Note that the reactions characterized by Lj„t are catalytic, whereas those of L/ ong and L s ^ ort 
are inhibitive. 

To incorporate the effects of these length scales and the inhomogeneities they produce, we introduce diffusion terms 
into the system of equations. Although diffusion is clearly appropriate for such things as the movement of material 
such as stars and clouds through the disk, it might seem more realistic to use a wave equation to represent the effects 
of Shockwaves and radiation. However, as is well-known, there are wavelike solutions to reaction-diffusion equations, 
such as the prototype Fisher equation pi] . Also, we can see from the success the cellular automata of Gerola, Seiden 
and Schulman || and the model of Elmegreen and Thomasson [^2| that this is a valid approach. 

Because the scale height of the Milky Way is much smaller than that of the radius of the disk, we consider the 
model on a two-dimensional annulus, leaving out the galactic bulge si nce i t has little influence on the star formation 



in the disk. Therefore the component functions we considered in Sec. II A are now functions of the radius p and the 
angle 9, in addition to time t: 
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c = c(p, 6; t) g — g(p,Q\t) s = s(p,9;t) r = r(p,9;t) h = h{p,6;t) 



Note that the radius and angle, along with the angular velocity u>(p) in the disk, will be the only Greek letters that 
are not parameters of the model. Being in a rotating system, with angular velocity we must use the convective 

derivative, so that the evolution of the components is given by 



dc 9c oiig 2 



UJ {P)'aa = a 2 Cr ~ Pich — 7 CS ~ e l c2 (2 a ) 



~H7 + ^yP) = \-a 2 cr + s + 70s + e 2 c +d (2b) 

ot oO r 

-^+Lu(p)-^=f3 2 ch + e 3 c 2 -s + D s V 2 s (2c) 

Qi +UJ (p)d0 =ms-Mc,9,r)+D r V 2 r (2d) 

— +uj(p)— = m s-cb 2 (c,9,h)+D s V 2 h (2e) 



where the parameters have the same meaning as before, in the one-zone model. 



III. NUMERICAL RESULTS AND CONCLUSIONS 



Now that we have our model, we want to know if there are any structures of the right size that will arise. To 
do this, we consider a linearization of the model, and see if there are any instabilities. We assume that the galaxy 
is two-dimensional to simplify the analysis, and expand equations ^ to linear order around the steady state (i.e. 
c = Co + C + • ■ ■ , and similarly for the other functions) . We assume instabilities with a wave vector k growing with time 
scale A, so that, for example, C = Cke xt cos(k • x). Note that we can find eigenfunctions of the differential operator 
made up of the Laplacian and the convective derivative . Then these equations give us a matrix equation of the 
form M^v = Xv a , where v a is a column vector made up of the component functions - v a = (C, S, G, R, H). We can 
solve this matrix equation for the eigenvalues A, and find which modes of instability are likely to grow exponentially 
as a function of the parameters. 

To carry out the linearization analysis, we pick an arbitrary set of parameters lying within the physical ranges, as 



described in Sec. II B. If we look at the maximum real eigenvalue (MRE), we get a sense of which unstable modes 
will grow at the fastest rate. In Figure 2, we see a graph of the MRE of the matrix M§ as a function of the logarithm 
of the wavenumber k — |k|, using the values 

a x = l\ = ryx/10 = IO772 = 1.0 a 2 = 0.005 f3 2 = 0.9 7 = 0.5 



S = 0.002 ei = 1.0 0i = 0.6 <j> 2 = 0.4 e 2 = 0.9 e 3 = 0.04 
and the three diffusion constants 

D h = 1 D s = 10~ 3 D r = 10 5 

Since we are using L- lnt (^ 10 2 parsecs) as our length scale, then we note several things about the value of the MRE. 
It is always negative - all modes are stable and will eventually decay to equilibrium, although some decay slower. 
These modes can possibly allow the formation of some spiral structure in the disk, even if it is not a permanent 
pattern. Note that the slowly decaying modes have pattern sizes of between 10 -2 Lj n t to 10 2 Li n t (1 to 10 4 parsecs). 
Unfortunately, this seems to be generic within the model, although there is the possibility that there exist parameter 
sets with positive MREs that we have not found. 

Once we had a set of parameters that, at least for some time, would form structures of the right size, the equations 
(||) were numerically simulated using finite differencing, with the time evolution given by operator splitting. Because 
the diffusion constant D r used in the linearized analysis is so much larger than and D s , it was decided to use 
the mean field approximation for the UV radiation: the radiation was spread instantaneously across the galactic disk 
between time steps, instead of diffusing. Because the time increment used in the simulation was on the order of the 
light crossing time of the disk (about 10 4 years), this is not too unphysical a proposition. To take into account the 
effects of gravity, a constant linear velocity v was given to all the material, approximating the situation in the Milky 
Way. This was done by using u = v/r in the convective derivative. 
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TABLE I. Description of the various components of the 
model; in the numerical simulation, these are all functions of 
radius, angle and time. 



Variable Description 

_ Cold gas in GMCs 

g Warm, ambient gas 

s Massive stars 

d Light stars 

r Density of UV radiation 

h Density of Shockwaves from supernovae 



We present a picture of our results in Figure 3. The initial data is given on an annulus with r m [ n = 50L; n t (5 x 10 3 
parsecs) and r max = 150L int (1.5 x 10 4 parsecs), and is just a gaussian "blip" for one component (w), at an arbitrarily 
choosen location. The disk is rotating with a constant linear speed v = 30L; n t/T, which is the velocity of the Sun in 
the Milky Way. The boundary conditions on the inner and outer edge of the annulus is such that the radial derivatives 
of all the components are zero, i.e. there is no flux through the edges. This comes from the rationale that there are 
no clouds or stars infalling from the outside, and the central bulge does not interact much with the annulus. The run 
covered a time of about 200r, a timespan of two billion years or approximately eight revolutions of the galaxy. As 
one can see from the figures, spiral structure indeed is formed. The spirals develop and grow in length as they rotate 
with the galactic disk. One could imagine with several such initial gaussian "blips" , several spiral arms would form, 
thus resembling a real, flocculant, spiral galaxy. 

However, after about lOOr or one billion years, the ends of a spiral arm, begin to interact (diffuse) with each 
other and eventually (at about 1.5 billion years) merge forming a ring like structure. This appears to be a generic 
limitation of the model. Several initial data and parameter sets were tried with the same results. It seems that only 
a modification of the model will be able to save the spirals from this eventual, unfortunate fate! Work on this issue is 
currently underway and shall be presented elsewhere. Even with this limitation, this model competes well with other 
current models. 

Beyond this issue, there are a number of future directions that this work can be taken. This run only used one 
specific set of initial data; one can hope that spiral formation does not depend strongly on initial conditions. We 
have made several runs with differing initial data, but we have not further explored this area. We are also lacking a 
detailed study of the effects of changes in the values of the parameters; because of the large number of parameters 
and components, however, it is difficult to see how an exhaustive search for a more fruitful choice could be made. 
Modifying the boundary conditions is another avenue to consider - for example, allowing radiation to flow in from 
outside the galaxy, or to have interaction between the disk and the central bulge. Also, we might alter the role of the 
UV radiation in the model. One consideration is that we have not given the radiation any type of radial fall-off as 
it leaves the source, but instead used a mean field approximation. This, unfortunately, heats the entire disk evenly. 
If, instead, there is a great deal of radiation in the space near the source, and little far away, this might give a more 
realistic feedback mechanism to inhibit star formation near the source but allow it further away. There are, of course, 
other ways that the model can be changed to make it more realistic, such as adding in more reaction mechanisms. 
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